On the efficient Monte Carlo implementation of path integrals 
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We demonstrate that the Levy-Ciesielski implementation of Lie- Trotter products enjoys several 
properties that make it extremely suitable for path-integral Monte Carlo simulations: fast computa- 
tion of paths, fast Monte Carlo sampling, and the ability to use different numbers of time slices for 
the different degrees of freedom, commensurate with the quantum effects. It is demonstrated that a 
Monte Carlo simulation for which particles or small groups of variables are updated in a sequential 
fashion has a statistical efficiency that is always comparable to or better than that of an all-particle 
or all- variable update sampler. The sequential sampler results in significant computational savings if 
updating a variable costs only a fraction of the cost for updating all variables simultaneously or if the 
variables are independent. In the Levy-Ciesielski representation, the path variables are grouped in a 
small number of layers, with the variables from the same layer being statistically independent. The 
superior performance of the fast sampling algorithm is shown to be a consequence of these observa- 
tions. Both mathematical arguments and numerical simulations are employed in order to quantify 
the computational advantages of the sequential sampler, the Levy-Ciesielski implementation of path 
integrals, and the fast sampling algorithm. 
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I. INTRODUCTION 

The availability of short-time approximations having 
fast asymptotic convergence 0, warrants a closer 
look to the Monte Carlo implementation of the resulting 
Lie- Trotter products. The superior convergence of path- 
integral methods employing such short-time approxima- 
tions is achieved under the assumption that the integra- 
tion against all path variables is performed in an exact 
fashion. In practical applications, this is never the case, 
except for low-dimensional problems. The efficiency of 
the methods suffers from the slow convergence of Monte 
Carlo integration. Indeed, if the convergence order of a 
certain technique is v, then the computational cost to 
achieve a given error e, as measured by the number of 
calls to the potential function, has the form 0] 



Costcx l/e 2+1 ^, 



(1) 



assuming that the Monte Carlo samples are independent. 
This formula demonstrates that we cannot defeat the 
slow convergence of the Monte Carlo simulation by in- 
definitely improving the convergence order. 

In practical applications, Eq. represents a very op- 
timistic evaluation, because one must deal with the ad- 
ditional problem of build-up of correlation among path 
variables, as the number of variables increases Q- Thus, 
only a small group of path variables can be updated in an 
efficient fashion at a time. Significant research has gone 
into the problem of diminishing the correlation between 
path variables and ensuring a more efficient sampling. 
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Techniques such as the staging method H, the thread- 
ing algorithm jaL the bisection method the multi- 
grid technique Jo], an d the normal mode and Fourier 
approaches 0, la, can significantly decrease the cor- 
relation times in path integral Monte Carlo simulations. 
A recent technique developed in Ref. IT?! which is called 
the fast sampling algorithm, builds upon some special 
properties of the so-called Levy-Ciesielski representation 
of the Feynman-Kac formula [llj . The technique can be 
regarded as part of the random series El (in the 
continuous form) or normal mode (in the discrete form) 
approaches to path integration. In the present work, we 
demonstrate that the technique is capable of reducing the 
computational time necessary to achieve a prescribed sta- 
tistical efficiency from n 2 calls to the potential function 
(scaling that is valid for most normal mode representa- 
tions) to n log 2 (n) . Here, n represents the number of 
path variables. 

In Section II, we analyze the computational cost of 
the Metropolis et al algorithm [lH Hij for high dimen- 
sional systems, from the point of view of statistical ef- 
ficiency. We present both mathematical and numerical 
arguments to justify the finding that updating particles 
one at a time is statistically at least as efficient as us- 
ing all-particle moves. The most important cases where 
updating particles or path variables one at a time re- 
sults in important computational savings are i) for clas- 
sical systems, the case where the computational time for 
the whole potential increases linearly with the compu- 
tational time necessary to update only one particle and 
ii) for path integral simulations, the case where the path 
variables can be grouped in independent random vectors. 

In Section III, we study the statistical efficiency of the 
Metropolis et al sampler for random series as well as for 
the normal mode implementation of Lie- Trotter prod- 
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ucts. We conclude that the Levy-Ciesielski representa- 
tion is superior in both cases, allowing for a reduction 
in the computational cost of log 2 (n)/n, by comparison 
with most normal mode implementations. The reader 
must realize that the fast sampling algorithm is not re- 
ally a sampling technique. Rather, it is a property of 
the Levy-Ciesielski representation and can only be uti- 
lized for path integrals. Similarly, the fast computation 
of paths is also a property of the Levy-Ciesielski series, 
rather than a technique. It enables the computation of 
paths in n\og 2 {n) operations instead of n 2 , the number 
necessary for most other normal mode implementations. 
It is true, starting from Coalson's Fourier-like normal 
mode approach Q , one can still construct the paths in a 
time proportional to n log 2 (n) , by using fast sine- Fourier 
transform. This has been observed by Mielke and Truh- 
lar [Tij . However, constructing a fast sampling algorithm 
by using this methodology is rather difficult. Perhaps the 
most important property of the Levy-Ciesielski series is 
that it constitutes a link between the continuous and the 
discrete path integral techniques ^lj. Thus, almost all 
algorithms developed for the discrete case have an ana- 
logue in the Levy-Ciesielski language. Li and Miller [l7| 
have recently demonstrated how a Lie- Trotter product 
for path integrals in many dimensions must be modified 
so that the number of time slices associated to each de- 
gree of freedom be proportional to the quantum effects. 
In the Appendix, we adapt the technique to the Levy- 
Ciesielski form and show how the number of time slices 
for each dimension should be chosen as a function of the 
particle masses. 

In Section IV, we utilize the fourth-order direct short- 
time approximation recently developed in Ref. |jj to ex- 
emplify the use of the Levy-Ciesielski series for the im- 
plementation of Lie- Trotter products. We then perform 
Monte Carlo simulations for the Nei9 Lennard-Jones 
cluster using both the all-variable update strategy and 
the fast sampling algorithm. The simulation is conducted 
at 4 K, using a number of 127 path variables per degree 
of freedom. The numerical results demonstrate that the 
standard deviations for the average energy and the heat 
capacity estimators are more than two times larger in 
the case when all particles are updated simultaneously. 
This translates in a computational saving of about 80% if 
the fast sampling algorithm is utilized. However, bigger 
computational savings are expected for larger numbers 
of path variables. 



II. CONSIDERATIONS ON THE STATISTICAL 
EFFICIENCY OF THE METROPOLIS SAMPLER 

In this section, we demonstrate that the maximal dis- 
placements in the Metropolis et al sampling algorithm 
decrease as fast as rt -1 / 2 with the number n of particles 
that are simultaneously updated. We argue that this de- 
crease is entropic in nature and has little to do with the 
interaction between particles. The generally accepted ex- 



planation for the decrease in the maximal displacements 
is that, by moving several of them at a time, we increase 
the chances that the particles collide. This explanation 
is mistaken and, to the contrary, we find that the de- 
crease in the acceptance probability exists even for non- 
interacting particles. By means of a numerical example 
we show that the entropic explanation also holds for par- 
ticles that interact through potentials having a strong 
repulsive part. 

Having quantified the decrease in statistical efficiency 
associated with multi-particle moves, we demonstrate 
that updating the particles one at a time (whether in a 
deterministic or random fashion) is the better strategy in 
terms of statistical efficiency. By statistical efficiency we 
understand the average distance covered by the random 
walker in the configuration space for a given computa- 
tional time and average acceptance probability. 

To begin with, let us assume that we are given a fi- 
nite collection X\, X 2 , . ■ . , X n of independent identically 
distributed random vectors (i.i.d.r.v's), taking values in 
some space M. d . These random vectors may represent, for 
instance, the space coordinates of a classical physical sys- 
tem made up of n identical particles that do not interact. 
Let p(x), with x S M. d , be the normalized distribution of 
any of the random vectors Xi. The distribution p(x) is 
assumed to be a smooth function, that is, to have contin- 
uous first order partial derivatives. Again, by referring to 
our physical system, if V(x) is the (common) potential 
in which the particles move, then we may set 



P(x) 



/QW), 



where /? is the inverse temperature and Q((3) is the con- 
figuration integral of the corresponding canonical system. 
By independence, the overall distribution of the ran- 
dom vectors is given by the product p(xi)p(x 2 ) . . . p(x„), 
which is a smooth distribution density on the space M. dn . 

It is perhaps clear that the best strategy for 
Monte Carlo sampling of the product distribution 
p(xi)p(x 2 ) . . . p(x„) is to perform the sampling individu- 
ally, for each random vector. Thus, following Metropolis 
et al [bH lllH , we propose a new position for the random 
vector Xi from the trial distribution T(yj|x,-), which is 
uniform in a d-dimensional hypercube centered about Xj 
and has maximal displacements A s , for s = 1,2, ... ,d 
(therefore the sides of the hypercube have lengths 2A S ). 
The move is then accepted with probability 



min < 1, 



p(y t )r(xj|y t ) 
p(x l )T(y l |x J ) 



(2) 



and rejected with the remaining probability. Repeating 
the procedure, one generates an ergodic Markov chain of 
stationary distribution p(x^). 

Undesirable high correlation between successive posi- 
tions in the Markov chain is the result of two factors i) 
high correlation in the proposal distribution, correlation 
that increases as the maximal displacements decrease, 
and ii) low acceptance probability. As a rule of thumb, 
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in order to minimize the correlation, one tunes the aver- 
age acceptance probability 



Ad 



x min < 1 



dx i dy i p(x i )T(y i \x, 
p(yO T ( x i|Yi) 



p( x T (yil x i 



(3) 



to a value of about 50%, by increasing or decreasing the 
maximal displacements, as appropriate [l5j |. 

Assume now that we sample the random vectors 
together and update all variables at once, using the 
trial distribution T(yi|xi) . . . T(y„|x„). The move to 
(yi ; y2, • ■ • ,Yn) is accepted with probability 



p(yQr(x,|yO 
p(x i )T(y i |x i ) 



(4) 



and rejected with the remaining probability. The average 
acceptance probability is given by the formula 



Ac n = I dxidyx--- dx n dy n p(xi)T(yi\xi) 

lR 2d 



■ p(xn)T(y n \x n ) min < 1, JJ 



p(y») T ( x »|y») 
p{x-i)T{yi\xi) 



(5) 



If one attempts such a strategy and utilizes the op- 
timal maximal displacements computed for the case of 
single particle moves, the average acceptance probability 
decreases according to the law [lCl ] 

-Hn 



Ac n 

where H > is the relative Shannon entropy 

>(y)T(x|y) 



(6) 



H = - 



s. 2d 



p(x)T(y\x) log 



p(x)T(y|x) 



dxdy. (7) 



To avoid such a catastrophic decrease in the accep- 
tance probability, we must decrease the maximal dis- 
placements, so that to minimize the Shannon entropy at a 
rate equal to l/n. More exactly, if H n is the Shannon en- 
tropy corresponding to new maximal displacements A SjW 
and Ac is the desired constant acceptance probability, 
then 



\og(Ac)/n. 



(8) 



We now show that, for sufficiently large n, the decrease 
in the maximal displacements is controlled by the Fisher 
entropy of the smooth distribution p(x). For a random 
variable (one-dimensional random vector), the estimate 
can be obtained as follows. We start with the approxi- 
mation 

dylog [p(x + y)/p(x)] = 



Hn 



dx 



p(x) 



dx 



'2Ai 



2Ai,„ 

•Ai 



dy log [l + p(x + y)/p(x)-l] (9) 



dx 



P( x ) 
'4Ai „ 



Ai 



dy[p(x + y)/p(x) - 1]" 



where we have retained the first non-vanishing term in 
the logarithm expansion. This approximation becomes 
exact in the limit of small Ai n . In fact, in the same 
limit, one may expand the density p(x + y) around the 
position x to first order and conclude that 



dx 



'4Ai, B 



Ai, 



Ai, 



dy [p{x + y)/p{x) - 1] 



dx 



p{x) 



Ai 



dy[p'(x)/p(x)] 2 y 2 (10) 
p' (x) 2 1 p(x)dx. 



4Ai,„ j_ A: 
A 2 



The last integral appearing in the preceding formula is 
recognized as the Fisher entropy of the smooth distribu- 
tion p{x). For cZ-dimensional spaces, by a similar argu- 
ment, the reader may obtain the general expression 



d 

-Va 

6 ^ 



2 

s . a 



[d s p(x)} /p(x)dx 



(11) 



By comparing Eq. Ijlll) with Eq. © , we conclude that 
the asymptotic scaling of the maximal displacements in 
the limit of a large number of particles or random vectors 
that are updated simultaneously is given be the formula 



A,. 



A" 



(12) 



Here, the quantities A° s are asymptotic constants that 
may have values slightly different from the optimal maxi- 
mal displacements A s for a one-particle or random- vector 
update. The decrease in the maximal displacements pre- 
dicted by Eq. I|12f> is somewhat unexpected, given that 
the particles do not interact. In fact, the usual explana- 
tion that the decrease in the maximal displacements for 
multi-particle updates is the result of an increased chance 
in collision does not hold under closer scrutiny. As for 
the case of independent particles, the decrease is solely 
an entropic effect. 

Rather than resorting to more sophisticated mathe- 
matics to demonstrate the entropic nature of the decrease 
in the maximal displacements, we give a numerical ex- 
ample, where we verify Eq. I112H by performing a Monte 
Carlo simulation in the classical canonical ensemble for 
the 19-particle Lennard-Jones cluster. We have employed 
the Neig implementation of LJ19. Although all Lennard- 
Jones clusters have essentially the same classical thermo- 
dynamics, as can be seen from employing reduced coor- 
dinates, we give here the exact parameters because, in 
the second part of the paper, we shall also use the Neig 
cluster for quantum simulations. 

The total potential energy of the Neig cluster is given 

by 



V, 



19 

E 

i<3 



19 



V LJ (r^)+E V c( r 



(13) 
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where Vlj (r^ ) is the Lennard- Jones potential describing 
the interaction between the particles i and j 



rij 



(14) 



and V c (i"i) is the confining potential 

- Ren 



V c (r ; ) = e LJ 



(15) 



The role of the confining potential is to prevent the evap- 
oration of the cluster, for the cluster by itself is not ther- 
modynamically stable. The cluster is confined to its cen- 
ter of mass R cm by a polynomial potential that increases 
abruptly beyond the confining radius of R c — 2.25<7lj- 
The values of the Lennard- Jones parameters <jlj and clj 
used are 2.749 A and 35.6 K, respectively 0. The mass 
of the Ne atom was set to mo = 20.0, the rounded atomic 
mass of the most abundant isotope. 

The simulation has been conducted for 8 intermedi- 
ate temperatures arranged in geometric progression be- 
tween T m in = 0.15eLj and T max = 0.35eLj. To reduce 
the equilibration times of the Metropolis samplers, the 8 
statistically independent parallel replicas have been in- 
volved in periodical exchanges of confi gura tions, accord- 
ing to the parallel tempering algorithm USUI- For each 
n = 1, 2, . . . , 19, the simulation has consisted of 25 blocks 
of one million n-particle moves. The n particles partici- 
pating in a single move have been randomly selected from 
the 19 existing particles. After each block, the maximal 
displacements are decreased or increased so that the ac- 
ceptance probability for the last block is 50%, to a sta- 
tistical accuracy of about 0.5%. 

The quantities ^/nAi yn for the replicas of lowest and 
largest temperatures are plotted in Fig. ^ One can see 
that the asymptotic scaling predicted by Eq. ifT^ is re- 
spected to a very good degree. Similar plots for the re- 
maining 6 intermediate replicas show the same excellent 
agreement between the theoretical prediction and the re- 
sults of the simulation. It is therefore quite clear that 
the decrease in the maximal displacements is solely an 
cntropic effect that has nothing to do with an increase in 
the chances of collision. After all, whether we move one 
particle at a time or all particles together, the average 
interactions that any group of particles suffer must be 
the same, at least for a well-equilibrated simulation. 

The penalty for the decrease in the maximal displace- 
ments is an increase in the correlation between successive 
Monte Carlo steps, increase that is due to the pronounced 
correlation in the proposal step. If one updates each par- 
ticle individually in a deterministic or random fashion, 
then, on average after n p — 19 moves, the position of 
each particle is sampled from a distribution that spans a 
distance proportional to A s . Because of the decrease in 
the maximal displacements, it also takes about n p Monte 
Carlo steps for the all-particle strategy to achieve a simi- 
lar statistical efficiency, that is, to guarantee that each of 
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FIG. 1: Scaling of the maximal displacements with the num- 
ber of particles that are simultaneously updated for the LJ19 
cluster, at two different temperatures. Only every other com- 
puted values are marked on the plot. The thin lines have been 
added to help guide the eye toward the asymptotic region. 



the particles have been sampled from a distribution that 
roughly spans the same distance A s . 

To demonstrate the last assertion, let us look at the 
distances that are spanned by the random walker along 
some arbitrary direction, after N Monte Carlo steps. Ne- 
glecting the corrections that appear because the moves 
are not accepted with probability one (these corrections 
do not change the overall scaling, as long as the accep- 
tance probability is kept constant; in addition, not ac- 
cepting the moves with probability one further reduces 
the distances spanned by the walker as well as the statis- 
tical efficiency of the sampler, making our final conclusion 
even stronger), the position of the random walker along 
the direction s is 

N 

X StN = x s + A s>np ^(2u k - 1), 
fc=i 

where the quantities u k are independent random vari- 
ables uniformly distributed on the interval [0, 1]. For suf- 
ficiently large N, the sum J2k(^ u k ~ 1) has a Gaussian 
distribution of variance N/3 centered about the origin, 
as follows from the central limit theorem. Therefore, the 
random variable X s ^ is a Gaussian centered about x s 
and of variance n N/3. The average distance relative 
to the starting point spanned by the random walker is 



\ Z \(2nAl np N/3)- 1/ \-^( 2 <^ 3 ) 



dz 



= ^2/{3n)A s . np N 1 ' 2 = A° s ^2N/(37rn p ). (16) 

From Eq. (|16f) . we see that it takes N ~ n p Monte Carlo 
steps for the all-particle strategy to achieve a statistical 
efficiency comparable to that of the one-particle strategy, 
also after n p Monte Carlo steps. 
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FIG. 2: Average distances spanned by the random walker 
after n v Monte Carlo steps as a function of the number of 
particles n that are updated simultaneously. The error bars 
are less than half the size of the plotting symbols. Only every 
other computed values are marked on the plot. 



In fact, the statistical efficiency remains roughly the 
same no matter how many particles we move simultane- 
ously. To illustrate this by a numerical example, we have 
evaluated the average distances spanned by the Monte 
Carlo walker after n p steps, while simultaneously updat- 
ing groups of n = 1, 2, . . . , 19 randomly chosen particles. 
If Xn denotes the position of the walker at time N, then 
the average distance is 



ous maximal displacements cannot be determined during 
the simulation and have to be fixed a priori. 

Let us conclude this section by mentioning that the 
decrease in the maximal displacement of n -1 / 2 is solely 
due to the larger number of particles that are updated 
simultaneously. If the strength of the correlation in- 
creases with the number of particles, then the decrease 
in the maximal displacement for the multiparticle update 
is n -1 / 2 times the decrease in the maximal displacement 
for the one-particle update. For instance, consider the 
problem of sampling the distribution 



exp 



1 (s - xx) 2 + (a?i - x 2 ) 2 + ... + (x n - xx) 2 

2 a 2 /n 



(18) 

where a 2 — h 2 (3/mo. This distribution is encountered in 
the construction of Lie- Trotter products. The decrease 
in the maximal displacement for the all-particle update 
strategy is n -1 / 2 x n -1 / 2 , with the first factor due to the 
larger number of particles and the second factor due to 
the decrease in the maximal displacement for one-particle 
moves. Thus, the number of Monte Carlo steps necessary 
to achieve a prescribed statistical efficiency is n 2 . Since 
the computational effort for a single Monte Carlo step in 
terms of calls to the potential is also proportional to n, 
we see that the cost for the direct Monte Carlo sampling 
of the Trotter-Lie products is proportional to n 3 , result 
consistent with the one obtained in Ref. 0- 



JV-l 



(\\X np - X \\) 



lim 

N-kx. 



N ^ 



\X. 



n p -\-k 



(17) 



k=0 



In collecting the averages, one must discard all differ- 
ences ||X n — Xk\\ for which a parallel tempering swap 
has occurred at any Monte Carlo step between k and 
k + n p . The average distances are shown in Fig. 0and 
are seen to closely mimic the behavior of the maximal dis- 
placements. Thus, the one-particle and the multiparticle 
updating strategies have essentially the same statistical 
efficiency. 

In practical applications, updating the particles one at 
a time is almost always the winning strategy in terms 
of computational effort for same statistical efficiency. In 
many applications, the potential can be decomposed in 
Up smaller parts, each describing the interaction of a 
particle with its environment, and each taking n p -times 
less computational effort to evaluate. Thus, the single- 
particle strategy ensures that each particle is sampled 
from a distribution spanning a distance of A s , in a time 
roughly equal to the time for a single all-particle update. 
The all-particle strategy takes n p -times more computa- 
tional resources to achieve similar results. But even for 
the situations where such a decomposition is not possi- 
ble, the sequential sampling is superior because it allows 
for a better tuning of the maximal displacements. In the 
all-particle strategy, the optimal ratios between the vari- 



III. STATISTICAL EFFICIENCY FOR PATH 
INTEGRAL SAMPLING 



In the preceding section, we have demonstrated that 
the efficiency of the sequential sampler cannot be de- 
feated by employing all-particle moves. This ffiiding sim- 
plifies the efficiency study for the different path-integral 
sampling strategies. In this section, we shall analyze 
the statistical efficiency of the random series approach to 
path integrals and of related implementations. We shall 
see that the computational time for a given statistical ef- 
ficiency scales as n times the cost to evaluate the action, 
for most series. Here, n is the number of path variables. 
One important exception is the Levy-Ciesielski series, for 
which the scaling is log 2 (n) times the cost to evaluate the 
action. This remarkable property of the Levy-Ciesielski 
series constitutes the engine behind the fast sampling al- 
gorithm [Toj | . 

In the second subsection, we specialize the findings ob- 
tained in the case of random series for the normal mode 
implementation of Lie- Trotter products. For such prod- 
ucts, the time to evaluate the action is proportional to 
the number of path variables. Therefore, the computa- 
tional time for a given statistical efficiency is proportional 
to n 2 for most normal mode approaches, except for the 
Levy-Ciesielski one, for which the scaling is n log 2 (n) . 
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A. Random series implementation of path integrals 

A standard approach for the numerical implementation 
of the Feynman-Kac formula 0, 12^, is via random 
series [lj, uM- The implementation is as follows. Let 
{Afc(r)}fc>o be any orthonormal basis in L 2 [0, 1] such that 
Aq(t) = 1. Define the primitives 



Afe(w) = / A fc (r)dr. 



(19) 



Let VI denote the set of all sequences a :— {ai,a 2 , . . .}. 
The Gaussian measure 



dP[a] = Yl 



1 



k=l 



2tt 



- al,2 da k 



(20) 



on Q makes the normal random variables a :— 
{0,1,02,... } independent and identically distributed. 
With the above notations, the one-dimensional Feynman- 
Kac formula reads 01 



p(x,x';f3) = p fp (x,x';P) / dP[ 



x exp I -0 J V 



c r (u) +gy^a fc Afc(u) 



k=l 



du } . (21) 



Eq. (I2H is called the random scries representation of 
the Feynman-Kac formula. The quantities p(x, x'; [3) and 
Pf P {x, x'\ (3) represent the density matrices of the physical 
system and of the free particle, respectively. x r (u) stands 
for x + {x 1 — x)u, whereas a = (h 2 /mo) 1 ' 2 ■ The gen- 
eralization to many dimensions is straightforward: one 
just considers an independent random series for each ad- 
ditional physical degree of freedom. 

The random series representation of the Feynman- 
Kac formula is made possible by the Ito-Nisio theorem 
[ill. |24| | . which gives an explicit construction of the Brow- 
nian bridge entering the Feynman-Kac formula. This 
theorem implies that the Feynman-Kac formula is in- 
variant to orthonormal transformations corresponding to 
changes from a basis {Xk(u)}k>i orthogonal on Aq = 1 



to another basis {Aj.(u)}fc>i, also orthogonal on the con- 
stant function. The reader may easily rationalize this 
observation by noticing that the measure defined by 
Eq. (|20|l is invariant under an orthonormal transforma- 
tion a' k = V, ,^.,-0. 

Important examples of series representations of the 
Feynman-Kac formula are provided by the Wiener- 
Fourier series and Levy-Ciesielski series. The Wiener- 
Fourier series representation is obtained from the cosine 
Fourier basis {Afc(-r) = \f2 cos(fc7rr)}fe>i, which, together 
with Ao(t) = 1, forms a complete orthonormal basis of 
L 2 [0, 1]. The primitives of the cosine functions are 



2 sm(kiru) 



A fe (u) = / \ k (T)dr 



JO _ V 7T' 

Upon replacement in Eq. J2TJ, we obtain 



p(x,x';(3) = Pf p (x, x';0) J dP[a] exp <j - f3 



x / V 
'0 



k=i 



2 sin(fc7ru) 



du}. (22) 



Eq. H22fl has been first utilized in the context of path in- 
tegrals by Doll and Freeman 12]. As an application that 
has no analogue for discrete path integral techniques, 
they have observed that the random series representation 
enables a computational technique, called partial averag- 
ing [25| , that has been recently shown to converge for all 
physically reasonable potentials 26]. 

A second important random series representation is the 
Levy-Ciesielski one. In this case, one starts with the so- 
called Haar basis, which is made up of the functions 

! 2 (fc-i)/2 ) r G ](Z — 1) /2 fe . I /2 k ] 
-2( fe -D/2, re [l/2 k ,(l + l)/2 k ] (23) 
0, elsewhere, 

where I = 2j — 1. Together with fo = 1, these func- 
tions make up a complete orthonormal basis in L 2 ( [0, 1]). 
Their primitives 



2 (fc-i)/2[ M _ (| _ l)/2% ue[(l- l)/2 fe , l/2 k ] 
Fk,j(u) = { 2( fe - 1 )/ 2 [(Z + l)/2 k -u], u G [l/2 k , (I + l)/2 k ] (24) 
0, elsewhere 



are called the Schauder functions. As McKean puts it 
|27j . the Schauder functions are "little tents," which can 
be obtained one from the other by dilatations and trans- 
lations. In modern terminology, this has to do with the 
fact that the original Haar wavelet basis is a multireso- 



lution analysis of £ 2 ([0, 1]) organized in "layers" indexed 
byfcj2|. 

For k = 1,2, . . . and j = 1,2,..., 2 k ~ 1 , the Schauder 
functions Fkj(u) are generated by translations and di- 
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for k > 1 and 1 < j < 2 



A— 1 



0.50 




FIG. 3: A plot of the renormalized Schauder functions for 
the layers k — 1, 2, and 3, showing the pyramidal structure. 



If we multiply them by 2~( fe ~ 1 )/ 2 j the Schauder func- 
tions make up a pyramidal structure organized in layers 
indexed by k, as shown in Fig. |3| The supports (the sets 
on which the functions do not vanish) of the Schauder 
functions are the open intervals of the form (ukj-i,Uk.j), 
for 1 < j < 2 fc_1 , where u^j = j2~( k ~ 1 \ The supports 
are disjoint for functions corresponding to the same layer 
k. Because of this property, we have the equality 



2 K-L 



— ak,[2 k ~ 1 u]+iFk,[2 k - 1 u]+i( u ), (27) 



latations of the function 



u, u€ (0,1/2], 
1-u, ue (1/2,1), 
0, elsewhere. 



More precisely, we have 

F k>j (u) = 2-( fe - 1 '/ 2 F 1 , 1 (2 fe - 1 W - j + I), 



(25) 



(26) 



for any sequence of numbers ak,i, a k ,2, ■ ■ ■ , flfc.2 fc - 1 • Here, 
[x] denotes the largest integer smaller or equal to 
x, whereas for u = 1, the quantities Ofc^-i+i an d 
-Pfc,2 fc - 1 +i(l) are defined to be equal to 0. 



In the new representation, the Feynman-Kac formula 
reads 0] 



p(x, x' ; 0) = Pf p {x, x 1 ; (3) I dP[a] exp ■ 

Jf2 



/ V 



x r {u) +o-^2a l[2 i-i u]+1 F l[2 i-i u]+1 (u) 
1=1 



du 



(28) 



r 



In the Levy-Ciesielski representation, the independent 
random variables 01,02,... have been re-indexed as 



{ai,j]l = 1,2, 



■J 



1, 2, . . . , 2 l 1 }, in agreement with 



the indexing scheme employed for the series. 

The numerical advantages of the Levy-Ciesielski rep- 
resentation are multiple. Assume that we truncate the 
series up to a number of n — 2 k — 1 path variables. That 
is, we use exactly k complete layers. Given u € [0,1], 
we only need k — \og 2 {n + 1) operations to perform the 
evaluation of the series at the point u. This is in con- 
trast with the Wiener- Fourier series, for which one needs 
n operations. This property is called fast computation of 
paths H2. 

A second property, which is called the fast sampling 
property 0, has to do with the sampling of the paths. 
We have already demonstrated in the preceding section 
that the efficiency of the sequential sampling technique 
cannot be exceeded by the techniques employing multi- 
variable updates. Also, notice that the maximal displace- 



ments for the individual update of the different path vari- 
ables, although not equal, do not decrease to zero. For 
path variables of large indexes, they converge to the max- 
imal displacements for a normally distributed random 
variable. (This observation is also true for the normal 
mode representation of Lie- Trotter products, considered 
in the following section). Since a complete sweep through 
the space of path variables is done in n steps, it follows 
that the computational effort to achieve a prescribed sta- 
tistical efficiency for the Wiener-Fourier series is n times 
the cost to evaluate the action (we shall call action the 
one-dimensional integral over the interval [0, 1] appear- 
ing at the exponent). This cost does not change if an 
all- variable sampling strategy is adopted. 

For the Levy-Ciesielski series, one still needs to update 
each path variable one at a time. However, it is not 
necessary to compute the whole action in order to update 
a variable. More precisely, if the variable a\j is to be 
updated, then one only needs to compute the quantity 
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r «2- (! - i) 

~ a '.j exp \ -p I V 

{ ./(j-x)2-«-i) 



a i,[2 1 - 1 M] + l-f),[2 ! - 1 «] + l(' lt ) 
Z=l 



(29) 



in order to make the decision if the variable aij is to be 
updated or not, according to the Metropolis et al crite- 
rion. The point here is that the terms 



exp 



(j-i)2-P-« 



V 



X r {u) 



k , . 

^ 0't,\2i-^u\+lPl,\2'- 1 u]+l (u) du> 
1=1 1 ' 



and 



exp 



(3 I V 



5>L 

(=1 



[2' 



av(w) + a 

] + l F L[2'- 1 u] + l( u ) 



(III 



do not contain the variable aij, because the functions 
Fij(u) vanish outside the open interval 



(( J --l)2-P- 1 ),j2-P- 1 )) 



Therefore, with a single evaluation of the action, we can 
update all 2 i_1 variables from the layer I, independently. 
Thus, the computational cost is the product between the 
number of layers k = log 2 (n + 1) and the cost to eval- 
uate the action. The Levy-Ciesielski representation is 
n/log 2 (n + 1) times faster than the Wiener- Fourier se- 
ries from the point of view of sampling efficiency 



B. Sampling efficiency for the normal mode 
approach to Lie- Trotter products 

The traditional way of constructing approximations 
to the Feynman-Kac formula is via Lie- Trotter products 
0, H3 . Such a construction starts with a short-time 
high-tcmperature approximation to the density matrix, 
say po(x, x'; (3). Because the density matrix of a free par- 
ticle is strictly positive, any such short-time approxima- 
tion can be put in the product form 



p (x, x'\ (3) = pf p (x, x'; (3)r (x, x'; (3). 



(30) 



Letting xq = x, = x' , and Ui — i/(n + l) for < i < 

n + 1, the n-th order Lie- Trotter product obtained from 



the short-time approximation considered above takes the 
form 

~ n 

p n (x,x';/3) = / Y\p a 2 (ui _ Ui+l) (xi,x i+ i) 



n 

x II r Q (x :j ,x : j + i;f3/2 k )dxi ■ ■ ■ dx n . 

3=0 



Here, p u {x,x') is defined by 

p u (x,x') = (27rw)- 1/2 exp [-(x' -xf/{2u)} 



(31) 



(32) 



A set of n Gaussian random variables having joint 
probability distribution 



Y[Pa 2 {u z -u z + 1 )(Xi,X i+ i)dxi ...d,X n 



(33) 



can be constructed in various ways WL 0, 0OJ . For in- 
stance, by diagonalization, Coalson 0] has shown that 
if <zi, 02, • • • , a n are independent Gaussian variables of 
mean zero and variances 



Ai = Aa 2 (n+ l)sin 2 
then 



2(n+l) 



1 < i < n, (34) 



X r (Ui) -)- ^ CLjSij 
3=1 

has the distribution given by Eq. I|33|) . Here, 

ijir 



(35) 



Si. 



sin 



, 1 < i, 3 < n. 



n+1 + 
As argued in R,efs. lTlIIT3l a more useful form is 



c r (ui) + a-y^ajSij/X 



1/2 



(36) 



(37) 



3=1 



with the variables ai, a%, . . . , a n being independent iden- 
tically distributed normal random variables. The main 
reason is that the temperature dependence is now buried 
into the expression of the short-time approximation. In 
the limit of large n, the resulting discrete approximation 
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Pn {x,x';P) = Pfp {x,x';p) f da,-- - [ da n {2^)- n/2 TT eT°> 1 

Jm jr k=1 



x II r » 

i=0 



-(ui) + a^S^jaj/X 1 / 2 ,x r (u i+1 ) -j-a^Si+ijaj/X] 



/a. 



n + 1 



(38) 



converges to the Feynman-Kac formula in re-scaled form. 
Therefore, the thermodynamic estimators obtained from 
formal differentiation against the inverse temperature 
have finite variance in the limit of large number of path 
variables [lji|^|22|]. 

As for random series, Eq. (|38|) is invariant under or- 
thogonal transformations. In the present form, the for- 
mula looks like the Wiener-Fourier series. However, as 
argued in Ref. 

in 

by appropriate orthogonal transforma- 
tions, the representation given by Eq. (|38|l can be made 
to look like any series we want [more precisely, in the 
limit of large n, we can make Eq. <|38[) look like any se- 
ries allowed by the Ito-Nisio theorem]. If n = 2 k — I, 



another possible construction of a set of n Gaussian vari- 
ables having the joint distribution given by Eq. (|33|) is 

k 

x r {uj) + cr^2a l[2 i-i u . ]+1 F lt[2 i-i u . ]+1 (uj). (39) 
1=1 

The exact orthogonal transformation that takes Eq. I|37|l 
into Eq. I|39|l does not really matter, as Eq. I|39|l can be 
demonstrated directly from the Le vy-C iesielski represen- 
tation of the Brownian bridge 0, [n| . The Lie- Trotter 
product now becomes 



p n (x,x';/3) = pfp(x,x';(3) / da^x ■ ■ ■ 

k 



k T 



da k>2k -i (2tt)-"/ 2 [] [] cxp (-al/2) 
i=i i=i 

r {u ) + a y^a;j2'- 1 M 1 1 + l^j2'- 1 t 1 ,1+l( M j)i 
1=1 

k 

(u 3+ i)+aJ2 a i,[2 ! - 1 « J+ i]+i^ 7 i,[2 ! - 1 « J+ i]+i( ,i j+i); /^/2 



3=0 



(40) 



Eq. 1)40(1 has the same numerical advantages over Coal- 
son's sine-Fourier form as the Levy-Ciesielski represen- 
tation has over the Wiener- Fourier series: fast compu- 
tation and sampling of paths. The analysis performed 
in the preceding section carries over here in a simple 
form. The computation and sampling of paths is done in 
(n + 1) log 2 (n + 1) operations for the Levy-Ciesielski rep- 
resentation because there are n+ 1 slices. More precisely, 
for sampling, (n+ 1) log 2 (n+ 1) represents the number of 
calls to tq[x, x'\ P/2 k ] in order to update all path variables 
sequentially. For Coalson's sine-Fourier form, one needs 
(n + l) 2 operations two perform the sampling (whether 



sequential or all-variable updates are attempted) in or- 
der to ensure a given statistical efficiency of the sampler. 
The computation of paths also take (n + l) 2 operations 
if implemented directly. However, the computation of 
paths can be done in (n + l)log 2 (n + 1) operations by 
means of the fast sine-Fourier transform, as pointed out 
by Mielke and Truhlar [l^. Most likely, the orthogonal 
transformation that takes Eq. (|37|l into Eq. I|39|l is the 
one that enables the sine-Fourier transform algorithm. 

In the Levy-Ciesielski representation, the quantity that 
must be used to test if the variable aij is updated or not 
is 
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;'2 



k — l + l_ 



e-<*/ 2 n r ° 

j=(i-l)2 fc - ! + 1 



r (Uj) + a y^a;,[2'-iu,l + iJ)j2'-iu,l + i(^j, 
1=1 

x r (uj+i) + cr^ ai,[2 ! -i« J+ i]+i-Fi,[2'-it lj+ i]+i( w j+i);/3/2 fc 



1=1 



(41) 



If the short-time approximation introduces additional 
path variables, then these path variables are to be sam- 
pled separately. They can be grouped into 2 k indepen- 
dent subsets, which can be individually tested for accep- 
tance. We shall give such an example in the following 
section. 

Let us address the issue of overall efficiency for the 
Monte Carlo simulation of Lie- Trotter products. Assume 
the order of convergence of the short-time approximation 
is v. To achieve a final error of e, we need to utilize 
n cx £~ x l v path variables. The computational cost to 
efficiently update all path variables once is proportional 
to n 2 cx £~~ 2 / v for the Wiener-Fourier approach and to 
nlog 2 (n) cx e -1 /"^ -1 log 2 (l/e) for the Levy-Ciesielski 
form. This cost must be multiplied by the number of 
steps necessary for the Monte Carlo sampler to reach an 
error of e, number of steps that is proportional to e~ 2 . 
Thus, the overall cost is 



Cost cx e- 2 - 2/ ", 
for the Wiener-Fourier approach, and 



Cost cx iz-V 3 - 1 /" log 2 (l/e), 



(42) 



(43) 



for the Levy-Ciesielski approach, respectively. The scal- 
ing for the Levy-Ciesielski approach is only marginally 
worse than the ideal scaling expressed by Eq. £[J. 



TABLE I: Quadrature points and weights for the 4-point 
Gauss-Legendre technique on the interval [0, 1]. 



i 


1 


2 


3 


4 


0, 


0.069431844 


0.330009478 


0.669990522 


0.930568156 




0.173927423 


0.326072577 


0.326072577 


0.173927423 



The fourth-order short-time approximation is given by 
the formula 



r (x,x';(3) = / (27r)-^ 2 e - (b i +b ^^ 2 

f 4 

x exp ^ -/3^2u k V 



k=l 



x r (9 k ) + aJ2 b j^j( e k) 



(44) 



In Eq. l|4"4"|) . uj k and 6 k are the weights and points for 
the four-point Gauss-Legendre quadrature technique on 
the interval [0, 1]. They are given in Table [Q for ease of 
reference. The three functions Aj(u) are defined by the 
equations 

Ai(u) — v3u(l — u), 

A 2 (u) = r(u)cos[ai(u- 0.5) + a 2 (u- 0.5) 3 ], (45) 
A 3 (u) = r(u) sin[ai(u - 0.5) + a 2 (u - 0.5) 3 ], 



IV. AN APPLICATION OF THE 
LEVY-CIESIELSKI IMPLEMENTATION 

In this section, we illustrate the numerical advantages 
of the Levy-Ciesielski implementation and the fast sam- 
pling algorithm by utilizing the technique in the context 
of the direct fourth-order short-time approximation in- 
troduced in Ref. 0. The resulting path-integral expres- 
sions are then employed to compute the energy and the 
heat capacity of the Neig cluster, at the temperature of 
4 K. We perform two simulations using the fast sam- 
pling algorithm and the all-variable sampling strategy. 
Important reductions in the statistical errors of the ther- 
modynamic energy and heat capacity estimators are ob- 
served for the fast sampling algorithm. These reductions 
are solely due to the decrease in correlation between the 
successive steps of the generated Monte Carlo Markov 
chain. 



with 

r(u) = {u(l -u)[l- 3u(l - u)]} 1/2 . 

The numerical values of the constants a.\ and a 2 are 

ai re 6.379716466 and a 2 w 8.160188248. (46) 

Using Eq. (|40|l . we can arrange the additional path 
variables as supplementary layers in the Levy-Ciesielski 
series. Extend the functions {A;(it);l < / < 3} outside 
the interval [0, 1] by setting them to zero and define 

G«(«) = 2- fc / 2 A«(2 fc tt-i + l) ) (47) 

for 1 < I < 3 and 1 < j < 2 fe . Then, with the convention 
that a/ ;2 '-i+i = 0, for 1 < I < k, and &j 2 fc +i = 0, for 
1 = 1,2, 3, we have Q 
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k T 



Pn ^ X '] P l = { da hl ...[ da k ^ (2.)-/ 2 op f - i 53 E a h 
p fp (x,x';/3) J R Jr \ 2 J~!~! 



x / dbi 



3 2" 

2 "i>3 
1=1 3'=1 



x exp — /3 y V 



{=1 

^XXpH+l G fe?[2*u]+l( U ) 



(«) 



(48) 



;=i 



The action integral is performed by means of the quadra- 
ture scheme specified by the 4 • 2 k — 4(n + 1) quadrature 
points 

u i:j = 2~ k {Q l +3-1), 1 < i < 4, 1 < j < 2 k (49) 
and the corresponding weights 

Wij = 2~ k Ui. (50) 
The quantities 9i and u>i are those from Table [I] 
I 



du 



r 



The additional path variables &;j make up three dif- 
ferent layers that are additional to the layers made up 
by the Lie- Trotter path variables aij. Such a layer I 
is selected randomly with probability equal to the other 
layers. Again, due to the fact that the functions Gij(u) 
vanish outside the interval ((j — l)2~~ k , j2~ k ), the vari- 
ables bij from a given layer I = 1,2,3 can and must be 
updated independently. The appropriate weight is given 
by the formula 



" b ^ /2 expi -fi^WsjV 



X r (u so ) +Cr5Z a J,[2 J - 1 u^]+l Fl\2 l - 1 u sj ] + l{ u sj) 
1 = 1 



1 = 1 

r 



For a Lie- Trotter variable a^j, the appropriate weight is 



(51) 



e ar .' /2 exp' 



j=l + (i~l)2 k - l + 1 s=l 



{=1 

3 

+0-XX[2* ! «„d+l G fc,[2*„ sJ ]+l^ 



(52) 



I 

Our choice of the Neig cluster for numerical experi- ber of path variables that is large enough to facilitate 

ments is motivated by the fact that the cluster presents the comparison, we conduct our computations at the low 

a deep classical global minimum [^H that is not destroyed temperature of 4 K. Numerical experiments show that, 

by the quantum effects j^. I n order to utilize a num- at this temperature, the number of path variables that 
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ensures a systematic error comparable to the statistical 
errors is 127 per degree of freedom (corresponding to the 
Trotter index n = 31). As opposed to the classical sim- 
ulation presented in Section II, here we do not use par- 
allel tempering to improve the sampling. This and other 
techniques that are commonly used to improve the qual- 
ity of the sampling have the property that they reduce 
the correlation of the Metropolis walker. Clearly, we do 
not want to measure the ability of the parallel tempering 
technique to do so. Therefore, we shall only conduct a 
simple Metropolis sampling, for each simulation. To en- 
sure that the basin associated with the global minimum 
is adequately sampled, we start all simulations from the 
global minimum. Had wc utilized a cluster with a dou- 
ble funnel topology of the potential electronic surface, 
we could not have ensured ergodicity of the simulation 
at 4 K. 

Our test consists in the evaluation of the average en- 
ergy and the heat capacity of the cluster using simultane- 
ous updates of all path variables and the fast sampling al- 
gorithm, respectively. The energy and the heat capacity 
estimators employed are those obtained by formal differ- 
entiation of Eq. 148( 1 . The estimators have been reviewed 
elsewhere [32J, |33| . In the case where all path variables 
are updated simultaneously, we have employed the same 
maximal displacements for all variables aij and bi t j. An 
optimal ratio between the maximal displacements for the 
physical coordinate x and the path variables has been 
determined in a separate Monte Carlo study, in which 
the two sets of variables were updated separately. It is 
however not at all clear if the ratio obtained this way 
remains the optimal one when all variables are updated 
simultaneously. This observation underlies again the bet- 
ter statistical efficiency of the sequential update: at the 
very least, the optimal displacements can be determined 
separately for each variable or group of variables during 
the original Monte Carlo simulation. Because for equi- 
librium properties the paths are closed (i.e., x' = x), the 
path variables from a same layer have identical marginal 
distributions and, therefore, they have identical maximal 
displacements. Thus, for the fast sampling strategy, a 
number of only 1 + log 2 (ra + 1) + 3 — 9 maximal displace- 
ments must be optimized. 

In both simulations, we have updated the coordinates 
associated with a given particle sequentially. That is, we 
randomly choose a particle and either update all vari- 
ables (for the all-particle update strategy) or only the 
variables associated with a randomly chosen layer (for 
the fast sampling strategy). In both cases, the computa- 
tional effort, as measured with respect to the number of 
calls to the potential, is the same. The simulations have 
consisted of 50 blocks of 20 thousand sweeps through the 
configuration space. Each simulation has been preceded 
by a number of 25 equilibration blocks. The statistical 
tests described in Ref. |32J have been employed to test for 
the independence of the block averages. 

The results of the two simulations are summarized in 
Table [H] The two sampling strategies have resulted in 



TABLE II: Energies, heat capacities, and associated statisti- 
cal errors (twice the standard deviation) for the Nei3 cluster, 
with the energy expressed in units of thj- 



Type of sampling 


Energy 


Heat capacity 


fast sampling 


-45.969 ± 0.023 


0.228 ± 0.069 


all variables 


-45.962 ± 0.052 


0.161 ±0.146 



similar values for the average energy. However, the sta- 
tistical errors are different, due to the higher correlation 
in the all-variable sampler. The energy statistical errors 
for the all- variable strategy are 2.26 times larger than 
those for the fast sampling strategy. It follows that the 
fast sampling algorithm allows for a saving of about 80% 
in the computational effort. For the heat capacity of the 
cluster, the saving is about 78%. For this particular ex- 
ample, the fast sampling algorithm has made the differ- 
ence between obtaining a reliable heat capacity and not. 
As discussed in the preceding section, for larger numbers 
of path variables, larger savings in the computational ef- 
fort are expected. Sure enough, the exact percentage de- 
pends not only on the quality of the sampling, but also 
on the smoothness of the estimator utilized. 



V. SUMMARY AND CONCLUSIONS 

The Monte Carlo sampler that uses all-particle or all- 
variable updates is not superior to the sequential sampler 
from the point of view of statistical efficiency. In fact, 
the sequential sampler is computationally more efficient 
whenever updating a single particle costs less than up- 
dating the whole potential. If one updates more than one 
particle at a time, the maximal displacements decrease 
inverse proportionally to the square root of the number 
of particles that are updated simultaneously. This effect 
has an entropic nature and appears even for independent 
variables. However, we warn the reader that the sta- 
tistical efficiency is only one factor controlling the rate 
of equilibration. The other factor is the rate at which 
the correlation in the generated Markov chain decays. 
As we have demonstrated in Section II, save the special 
case mentioned above, the all-particle update and the 
sequential samplers share roughly the same efficiency in 
terms of total volume in the configuration space that is 
sampled for a given computational effort. Which of the 
two techniques have a faster equilibration time when the 
statistical efficiency is the same has not been decided. 
Nevertheless, for independent variables, the sequential 
update is always superior. 

We have observed that, in the Levy-Ciesielski form, the 
path variables generated by the Lie- Trotter products are 
grouped in a small number of layers, with the variables 
from the same layer being statistically independent. This 
property, together with the observation that a sequential 
sampler has better statistical efficiency for independent 
variables, explains the superiority of the Levy-Ciesielski 
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representation versus the other random series or normal 
mode approaches. In the Levy-Ciesielski representation, 
by using a sequential sampler, one can efficiently update 
all variables using nlog 2 (n) calls to the potential. For 
most other normal mode approaches, the scaling is n 2 , 
whether a sequential or all-particle update sampler is uti- 
lized. 

To summarize, the computationally advantageous fea- 
tures of the Levy-Ciesielski implementation of path inte- 
grals are: fast computation of paths, fast path sampling, 
and the ability to use different numbers of path variables 
for the different degrees of freedom, commensurate with 
the quantum effects. The last property is discussed in 
the Appendix, where a relation between the number of 
time slices and the particle masses is suggested. These 
features recommend the Levy-Ciesielski representation as 
a useful technique for path integral implementations. 
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APPENDIX A: DIFFERENT NUMBERS OF 
TIME SLICES FOR DIFFERENT DEGREES OF 
FREEDOM 

Another advantage of the Levy-Ciesielski approach is 
that it allows for the utilization of different numbers of 
time slices for different degrees of freedom. Li and Miller 
|17| have recently shown how this must be done for gen- 
eral Lie- Trotter products. For the Levy-Ciesielski form 
of the fourth order short-time approximation, the Li and 
Miller procedure is equivalent to utilizing a larger num- 
ber of levels and quadrature points when particles with 
lighter masses are sampled. 

For definiteness, let us assume that the number of lay- 
ers are k and k', with k > k'. For the "light" coordinate, 
we utilize an entire random sum 

k 

x sj = x r (u sj ) + a 2J ai\2 l - 1 u s] ]+i F L[2 l - 1 u s] ]+i{ u s]) 
1=1 

3 

■^E^i^i+i G fc,W, j]+ iK)- ( Ai ) 
i=i 



where Uj — j/2 k , for j = 0, 1, . . . , 2 k . For the "heavy" 
coordinate, we utilize the independent sum 

ft' 

y S j = y r (u' sj ) +°"'XX[2<- 1 uy+i ^.p'-^y+iK-i) 
1=1 S1 

+^E 6 ;,[2^y +1 ^ [2 ^M +1 K J -) ! (A2) 
z=i 81 

where 

= 2^ k '(6 l + [j2 k '- k ] - 1), 1 < i < 4, 1 < j < 2 k . 

(A3) 

For example, if k — k' — 1, then the consecutive values 
u' s j and u' s J+1 for even j are equal. When computing the 
second term in the sum 

w sj V (x 3 j,y sj ) + w sJ+ iV (x 3 ,j + i,y St j + i) 

= w sj V (x sj ,y sj ) + w sJ+1 V (x S! j + i,y sj ) , (A4) 

one exploits the fact that, at least for many empirical 
potentials, it is easier to compute the difference 

V (x s>j+1 ,y sj ) -V(x sj ,y sj ) . (A5) 

Assuming that the strength of the interactions felt by 
each of the particles is the same, the number of levels that 
is appropriate for each degree of freedom is determined 
from the condition that the physical distances spanned 
by the variables from the last layers be equal. As can be 
seen from Eqs. (|A1|) . I|A2(1 . and H47|) . these distances are 
proportional to o2~ k l 2 and er'2~ fc / 2 , respectively. Re- 
membering that a = (h 2 (3/nio) 1 / 2 and a 1 = (^/J/mg) 1 / 2 
it follows that the numbers of layers k and k' must satisfy 
the relation 

m Q 2 k w m' 2 k ' (A6) 

or, equivalently, the number of slices n + 1 = 2 k and 
n' + 1 = 2 k must be in the relation 

(n + l)mo » (n' + l)m' . (A7) 

Eq. (|A7I) tells us that the appropriate number of time 
slices is proportional to the inverse mass of the particle. 
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